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Abstract 

The core motivating ideas of the space -time CESE method are clearly presented and critically analyzed. It is 
explained why these ideas result in all the simplifying and enabling features of the CESE method. A thorough 
discussion of the a scheme, a two-level non-dissipative CESE solver of a simple advection equation with two 
independent mesh variables and two equations per mesh point is also presented. It is shown that the scheme 
possesses some rather intriguing properties such as: (i) its two independent mesh variables separately satisfy two 
decoupled three-level leapfrog schemes and (ii) it shares with the leapfrog scheme the same amplification factors, 
even though the a scheme and the leapfrog scheme have completely different origins and structures. It is also 
explained why the leapfrog scheme is not as robust as the a scheme. The amplification factors/matrices of several 
non-dissipative schemes are carefully studied and the key properties that contribute to their non-dissipatedness are 
clearly spelled out. Finally we define and establish space-time inversion (STI) invariance for several non-dissipative 
schemes and show that their non-dissipatedness is a result of their STI invariance. 


1. Introduction 

The space-time conservation element and solution element (CESE) method is a high-resolution and 
genuinely multidimensional method for solving conservation laws [1-71]. Its nontraditional features include: 

(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a novel time marching strategy 
that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated 
without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability 
to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a 
non-dissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v) 
the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as 
independent marching variables to be solve for simultaneously. Note that CEs are non-overlapping space- 
time subdomains introduced such that (i) the computational domain can be filled by these subdomains; and 

(ii) flux conservation can be enforced over each of them and also over the union of any combination of them. 
On the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can 
be divided into several component parts with each of them belonging to a unique SE; and (ii) within a SE, 
any physical flux vector is approximated using simple smooth functions. In general, a CE does not coincide 
with a SE. 

Without using flux-splitting or other special techniques, since its inception in 1991 [1] the unstructured- 
mesh compatible CESE method has been used to obtain numerous accurate ID, 2D and 3D steady and 
unsteady flow solutions with Mach numbers ranging from 0.0028 to 10 [51]. The physical phenomena 
modeled include traveling and interacting shocks, acoustic waves, vortex shedding, viscous flows, detonation 
waves, cavitation, flows in fluid film bearings, heat conduction with melting and/or freezing, electrodynamics, 
MHD vortex, hydraulic jump, crystal growth, and chromatographic problems [3-71]. In particular, the rather 
unique capability of the CESE method to resolve both strong shocks and small disturbances (e.g., acoustic 
waves) simultaneously [13,15,16] makes it an effective tool for attacking computational aeroacoustics (CAA) 
problems. Note that the fact that second-order CESE schemes can solve CAA problems accurately is an 


NASA/TM— 2006-2 14385 


1 



exception to the commonly-held belief that a second-order scheme is not adequate for solving CAA problems. 
Also note that, while numerical dissipation is needed for shock capturing, it may also result in annihilation of 
small disturbances. Thus a solver that can handle both strong shocks and small disturbances simultaneously 
must be able to overcome this difficulty. 

In spite of its nontraditional features and potent capabilities, the core ideas of the CESE method are 
simple. In fact, all of its key features are the inescapable results of an honest pursuit driven by these 
simple ideas. The first and foremost is the belief that the method must be solid in physics. As such, in 
the CESE development, conservation laws are enforced locally and globally in their natural space-time unity 
forms for ID, 2D and 3D cases. Moreover, because direct physical interaction generally occurs only among 
the immediate neighbors, use of the simplest stencil also becomes a CESE requirement. Obviously, this 
requirement is also very helpful in simplifying boundary-condition implementation. 

The second idea is derived from the realization that stability and accuracy are two competing issues 
in time-accurate computations, i.e., too much numerical dissipation would degrade accuracy while too little 
of it will cause instability. In other words, to meet both accuracy and stability requirements, computation 
must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a real 
dilemma in numerical method development. As an example, schemes with high-orcler accuracy generally has 
high accuracy and low numerical dissipation. However, it is susceptible to instability. In fact, in dealing with 
complicated real-world problems, stability of these schemes often is difficult to maintain without resorting 
to ad hoc treatments. To confront this issue head-on, in CESE development, it is required that a solver 
be built from a non-dissipative (i.e., neutrally stable) core scheme. By definition, computations involving 
a neutrally stable scheme are performed right on the edge of instability and therefore the numerical results 
generated are non-dissipative. As such numerical dissipation can be controlled effectively if the deviation of 
a solver from its non-dissipative core scheme can be adjusted using some built-in parameters. Note that the 
above idea also plays an essential role in the recent successful development of a family of Courant number 
insensitive schemes [59,61,64,65,67]. 

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner; 
(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of 
ID schemes; (iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the 
method is compatible to the simplest unstructured meshes and thus can be used to solve problems with 
complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible, 
and special techniques that has only limited applicability and may cause undesirable side effects be avoided. 
Fortunately for the CESE development, as it turns out, the realization of the above lesser ideas (i)-(iv) 
follows effortlessly from that of the first two core ideas. 

From the above review, it is seen that the most critical task in the CESE development is to develop 
conservation-law enforcing and neutrally stable schemes. To give a preview over how this was accomplished, 
consider the simple advection equation 


du du 

m +a d^ = 0 


(i.i) 


where a ^ 0 is a constant. Let Xi = x, and X 2 = the considered as the coordinates of a two-dimensional Eu- 
clidean space E 2 . Then, because Eq. (1.1) can be expressed as V • h = 0 with h = ( au,u ), Gauss’ divergence 
theorem in the space-time E 2 implies that Eq. (1.1) is the differential form of the integral conservation law 


< j) h ■ ds = 0 (1.2) 

JS(V ) 

As depicted in Fig. 1, here (i) S(V) is the boundary of an arbitrary space-time region V in E 2l and (ii) 
ds = da n with da and n, respectively, being the area and the unit outward normal of a surface element 
on -S'(V'). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through the surface 
element ds, Eq. (1.2) simply states that the total space-time flux of h leaving V through S(V) vanishes; 
(ii) in E 2 , da is the length of a line segment on the simple closed curve 5(E); and (iii) all mathematical 
operations can be carried out as though E 2 were an ordinary two-dimensional Euclidean space. 
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It is well known that a solution to Eq. (1.1) (or Eq. (1.2)) represents non-dissipative data propagation 
along its characteristic lines defined by dx/dt = a. Moreover, as will be further discussed in Sec. 4, Eq. (1.1) 
is invariant under space-time inversion (STI), i.e., it transforms back to itself if x and t are replaced by —x 
and —t, respectively. Note that, in physics, STI invariance generally is referred to as PT invariance where P 
denotes parity (i.e., a mirror-image or spatial-reflection operation) while T denotes a time-reversal operation. 

From the above observations, it is obvious to the current author that non-dissipative behavior of the 
physical solutions to Eq. (1.1) will be preserved by the numerical solutions generated by a solver of Eq. (1.1) 
only if the solver also possesses the same STI invariance property. In other words, a non-dissipative solver of 
Eq. (1.1) must be STI invariant in some yet-to-be-defined sense. This conviction was further enhanced by the 
observation that the configurations of the stencils of all classical non-dissipative solvers of Eq. (1.1) such as 
the Wendroff [p.503, 72], the Crank-Nicolson [p.504, 72] and the leapfrog [p.100, 73] schemes do not change 
under space-time inversion while those of dissipative solvers of Eq. (1.1) such as the Lax scheme [p.97, 73] 
do change under space-time inversion. The preliminary results of a study along this line of thoughts were 
presented in a conference that took place in 1992 [2]. 

In a 1993 AIAA paper [74], the concept of STI invariance was also used by Thomas and Roe* to construct 
non-dissipative numerical schemes, (note: for the case where space-time is a 2D Euclidean space Ei, STI 
invariance is equivalent to the rotational symmetry referred to in [74]. However, in a Euclidean space of 
odd dimension, STI invariance and rotational symmetry are different because, in En, the determinant of the 
transformation matrix for STI is (— l) w while that for a proper rotation is always 1). In a 1994 ICASE report 
[75], the concept of STI invariance is again used by Roe to construct non-dissipative linear Bicharacteristic 
schemes. 

In addition to the non-dissipative property referred to above, a solution to Eq. (1.1) also possesses the 
following properties: (i) it is completely determined by the data specified at the initial time level; (ii) its 
value at a spatial point at a later time has a finite domain of dependence (a point) at an earlier time. As 
such, in the initial CESE development, the focus is on the construction of an ideal core solver of Eq. (1.1) 
that is flux-conserving and also possesses the properties similar to those listed above, i.e., a flux-conserving 
two-level, explicit and non-dissipative solver. For a reason to be presented in Sec. 4, none of classical solvers 
meets the above requirement. Fortunately, by utilizing the concept of STI invariance, we are able to construct 
such an ideal scheme which is referred to as the a scheme. It is a two level scheme with two independent 
mesh variables and two equations per mesh point. 

The rest of the paper is organized as follows: In sec. 2, a thorough discussion of the a scheme and its many 
equivalent forms is given. In addition, we establish a rather surprising result, i.e., the two independent mesh 
variables ul] and (u$)" of the a scheme separately satisfy two decoupled leapfrog schemes Eqs. (2.23) and 
(2.24) even though the a scheme and the leapfrog scheme have completely different origins and structures. We 
also explain why the leapfrog scheme is not as robust as the a scheme. In Sec. 3, we study the amplification 
factors/matrices of several non-dissipative schemes and point out the key properties that contribute to their 
non-dissipatedness. In Sec. 4, we define and establish STI invariance for each of the schemes discussed in 
Sec. 3 and show that their non-dissipatedness is a result of their STI invariance. Conclusions and discussions 
are given in Sec. 5. Finally, a computer code which can be used to verify Eqs. (2.23) and (2.24) numerically 
is listed in Appendix. 

2. The a scheme and the leapfrog scheme 

To proceed, consider the set fii of space-time staggered mesh points ( j,n ) (dots in Fig. 2(a)), where 

Hpf 

fii = {(j, ri)\j, n = 0, ±1, ±2, ±3, . . . , and (j + n) is an odd integer} (2.1) 

Each (j,n) G fii is associated with a solution element, i.e., SE(j, n). By definition, SE(j, n) is the interior of 
the space-time region bounded by a dashed curve depicted in Fig. 2(b). It includes a horizontal line segment, 
a vertical line segment, and their immediate neighborhood. 

* A participant of the 23rd Conference on Modeling and Simulation, April 30-May 1, 1992, Pittsburgh, PA. 
He was in the audience when the paper [2] was presented by this author, and participated in the discussion 
afterward. 
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At this juncture, the reader is warned that the notation used here may differ from that used in the 
previous CESE papers [1-71]. In particular, (i) the mesh indices j and n are only allowed to be whole 
integers here; and (ii) the spatial and temporal intervals that are denoted by ax/2 and At/2, respectively, 
in [1-71] are represented by Ax and At, respectively, here. These changes are introduced so that the CESE 
schemes and traditional schemes to be discussed here can be compared on the same footing. 

Let (x,t) € SE (j,n). Then Eq. (1.2) will be simulated numerically assuming that u(x,t) and h(x,t), 
respectively, are approximated by 

u*(x,t;j,n ) d = Uj + (u x )j(x - xj) + (« t )" (f - t n ) (2.2) 

and 

h*(x,t;j,n) = f (au*(x, t ; j, n), u*(x,t;j,n)) (2.3) 

Note that: (i) (u x )", and (w t )” are constants in SE(j,n), and they can be considered as the numerical 

analogies of the values of u, du/dx, and du/dt at the mesh point (j, n), respectively, (ii) ( Xj,t n ) are the 
coordinates of the mesh point (j, n) with Xj = j Ax and t n = nAt, and (iii) Eq. (2.3) is the numerical analogy 
of the definition h = ( au,u ). 

Let u = u*(x,t;j,n ) satisfy Eq. (1.1) within SE(j, n). Then one has (u t )" = — a(w x )". As a result, 
Eq. (2.2) reduces to 

u*(x,t;j, n) = u// + (u x )j [(a; - Xj) -a(t- t n )\ (2.4) 

i.e., u'j and (tt x )” are the only independent mesh variables associated with ( j,n ). 

Let E 2 be divided into non-overlapping rectangular regions (see Fig. 2(a)) referred to as basic conser- 
vation elements (BCEs). As depicted in Figs. 2(c)-2(e), (i) each (j, n) £ is assigned with two BCEs, 
i.e., CE_(j, n) and CE + (j, n); (ii) each BCE has one and only one pair of diagonally opposite vertices which 
belong to Oi; (iii) the space-time E 2 can be filled by CE_(j, n) and CE + (j, n), (j, n) £ fii; and (iv) CE(j, n), 
which is the union of CE_(j, n) and CE_|_(_ 7 , n), is referred to as a compounded conservation element (CCE). 
Let the flux of h* conserve over all BCEs, i.e., 


and 


h* ■ ds = 0, 


S(CE + U,n )) 


{j,n) £ fir 


(2.5) 


S(CE.{j,n)) 


h* ■ ds = 0, 


(j, n) £ Qi 


( 2 . 6 ) 


Note that, among the line segments forming the boundary of CE_(j, n), AB and AD belong to SE(j, n), 
while CB and CD belong to SE(j — 1, n— 1). Similarly, the boundary of CE + (j, n) belongs to either SE(j, n) 
or SE(j + l,n — 1). With the aid of the above observations and Eqs. (2.3) and (2.4), it can be shown that 
Eqs. (2.5) and (2.6) are equivalent to 


and 


(1 - v) [u + (1 + e)Ux}] = (1 - v) [u - (1 + v)u x \^ + l , 

{j,n) £ fii 

(2.7) 

(1 + v) [u - (1 - v)u s ]™ = (1 + u) [u + (1 - v)u x ]™zl I 

(j,n) £ Oi 

(2.8) 

(i) v = f a At/ Ax is the Courant number; (ii) 



(%)” d = 


(2.9) 


and (iii) to simplify notation, in the above and hereafter we adopt a convention that can be explained using 
an expression on the left side of Eq. (2.7) as an example, i.e., 


[it + (1 + v)Ux\j — «” + (1 + u){UxYj 
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At this juncture, note that: 

(a) Because 

du ax du 

dx 2 dx 

def 

if x = x/(ax/2), the normalized parameter (tij)" may be interpreted as the numerical analogy of the 
value at (j, n) of the derivative of u with respect to the normalized coordinate x. 

(b) Derivation of Eqs. (2.7) and (2.8) can be facilitated by the following observations: because u*(x,t;j,n ) 
is linear in x and t, it can be shown that the total flux of h* leaving CE _(j, n) or CE + (j, n) through 
any of the four line segments that form its boundary is equal to the scalar product of the vector h* 
evaluated at the midpoint of the line segment and the “surface” vector (i.e., the unit outward normal 
multiplied by the length) of the line segment. 

(c) By definition, point B depicted in Fig. 2(c) does not belong to either SE(j, n) or SE(j — 1 ,n — 1), i.e., 
h* is not defined there. However, this is not a problem for flux evaluation over AB or CB because the 
value of h* at an isolated point does not contribute to the flux of h* over a finite line segment. 

Let 1 — v 2 ^ 0, i.e. 1 — v ^ 0 and 1 + v ^ 0. Then Eqs. (2.7) and (2.8) reduce to 


U + (1 + v)u s } n 3 = [u - (1 + v)u s ] n j+ l , 

(j,n) G Hi 

(2.10) 

u - (1 - u)Ux]j = [u + (1 - X/)«x]"ri , 

(j,n) G Hi 

(2.11) 


respectively. Note that, for each (j, n) G Hi, Eq. (2.10) is a stronger condition than Eq. (2.7) (and therefore 
its equivalent Eq. (2.5)) in the sense that the former implies the latter for any given v while the latter implies 
the former only if an extra condition (i.e., v ^ 1 for this case) is imposed. Similarly, for each (J, n) G Hi, 
Eq. (2.11) is also a stronger condition than Eq. (2.8) and its equivalent Eq. (2.6). Hereafter, for any given 
i/, the system of equations defined by Eqs. (2.10) and (2.11) for all ( j , n ) G Hi is referred to as the first basic 
form of the a scheme. According to the above observation, any solution to this form automatically satisfies 
the conservation conditions represented by Eqs. (2.5) and (2.6). 

For any given (j, n) G Hi, Eqs. (2.10) and (2.11) form a pair of linear equations for the unknowns w™ 
and (%)" if the expressions on the right sides are treated as known source terms. It can be shown that, for 
any given v, the unique solution to this system is given by 

w" = ^{(i - ^) [u - (1 + ^)us]"+i + (1 + v) [u + (1 - zz)Mx]"ri}, ( j , n) G Hi (2.12) 

and 

{Ux)] = ^{[u- (! + ") u x]j+i - l u + (i - (j,n)Gfi 1 (2.13) 

In other words, for any given v and any (j,n) G Hi, Eqs. (2.10) and (2.11) imply Eqs. (2.12) and (2.13), and 

vice versa. Because u” and (u s )” can be explicitly evaluated in terms of * and (m®)"^ 1 using Eqs. (2.12) 
and (2.13), hereafter, for any given v, the system of equations defined by these equations for all ( j,n ) G Hi 
is referred to as the forward marching form of the a scheme. Obviously a solution to the forward marching 
form must also be a solution to the first basic form and vice versa. 

Next, by (i) replacing the indices j and n in Eq. (2.10) with j — 1 and n + 1, respectively, and (ii) using 
the fact that (j, n) G Hi <G> (j + 1, n — 1) G Hi (j — 1, n + 1) G Hi (hereafter the symbol <t=> is used as 

the shorthand for “if and only if”), one can show that, for any given v, the system of equations defined by 

Eq. (2.10) for all (j,ri) G Hi is identical to that defined by 

[u - (1 + zz)ux]" = [«+(! + i')u 5 ;]”j( 1 1 , (j, n) G Hi (2-14) 
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for all (j, n) G Hi. Similarly by (i) replacing the indices j and n in Eq. (2.11) with j+ 1 and n+1, respectively, 
and (ii) using the fact that (j, n ) G Hi <t=> (j — 1, n — 1) G Hi <G> ( j + 1, n + 1) € fij-, one can show that, for 
any given u, the system of equations defined by Eq. (2.11) for all (j, n) € Hi is identical to that defined by 

[u + (1 - v)ux\™ = [u - (1 - , (j, n) G Hi (2.15) 

for all (j,n) G Hi. Hereafter, for any given v, the system of equations defined by Eqs. (2.14) and (2.15) for 
all (j, n) G Hi is referred to as the second basic form of the a scheme. Obviously the first and the second 
basic forms are simply different representations of the same system of equations. 

For any given (j, n) G Hi, Eqs. (2.14) and (2.15) form a pair of linear equations for the unknowns tt” 
and (%)" if the expressions on the right sides are treated as known source terms. It can be shown that, for 
any given is, the unique solution to this system is given by 

Uj = ^{(1 + v) [u - (1 - i/)Mx ]"+ 1 + (1 - v) [u+ (1 + (j,n) G Hi (2.16) 

and 

(«*)" = [u - (1 - ^)ux]"+i - [u + (1 + (j, n) G Hi (2.17) 

In other words, for any given is and any (j,n) G Hi, Eqs. (2.14) and (2.15) imply Eqs. (2.16) and (2.17), and 
vice versa. Because it" and (%)" can be explicitly evaluated in terms of u"±i and (u $ )"£ using Eqs. (2.16) 
and (2.17), hereafter, for any given is, the system of equations defined by these equations for all ( j,n ) G Hi 
will be referred to as the backward marching form of the a scheme. Obviously a solution to the backward 
marching form must also be a solution to the second basic form and vice versa. 

The above developments are summarized and further elucidated in the following remarks: 

(a) According to the above discussions, among the four forms of the a scheme defined above, any one implies 
each of other three and vice versa. Because of this equivalency, hereafter the same essential conditions 
represented by any one of the above four forms will be referred to simply as the a scheme. 

(b) Because any solution to the first basic form of the a scheme automatically satisfies the conservation 
conditions Eqs. (2.5) and (2.6), it follows from the remarks made in item (a) that any solution to any 
other form also satisfies Eqs. (2.5) and (2.6). 

(c) Each component equation of the first or second basic form represents an implicit relation involving the 
mesh variables at two diagonally opposite neighboring mesh points, i.e., each “basic” stencil of the a 
scheme is formed by two neighboring mesh points. On the other hand, each component equation of the 
forward or backward marching form represents an explicit relation involving the mesh variables at three 
neighboring mesh points, i.e., each “marching” stencil of the a scheme is formed by three neighboring 
mesh points. 

(d) Because (i) the vector h* at any surface element lying on any interface separating two neighboring BCEs 
is evaluated using the information from a single SE, and (ii) the unit outward normal vector on the 
surface element pointing outward from one of these two neighboring BCEs is exactly the negative of that 
pointing outward from another BCE, one concludes that the flux leaving one of these BCEs through the 
interface is the negative of that leaving another BCE through the same interface. Due to this interface 
flux cancelation, the local conservation relations Eqs. (2.5) and (2.6) lead to a global flux conservation 
relation, i.e., the total flux of h* leaving the boundary of any space-time region that is the union of any 
combination of BCEs vanishes. In particular, because CE(j, n ) is the union of CE_(j, n) and CE + (j, n), 

f h*-ds = 0, (j,n)GHi (2.18) 

Js(CE(j,n )) 

must follow from Eqs. (2.5) and (2.6). In fact, it can be shown that Eq. (2.18) is equivalent to Eq. (2.12). 
Also it can be shown that the conservation condition over the union of CE + (j — 1 ,n + 1) and CE_(j + 
1 ,n+ 1) is equivalent to Eq. (2.16) (see Figs. 2(f)-(h)). 
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(e) As shown in [5], the two amplification factors of the a scheme are identical to those of the leapfrog 
scheme. As a result, the a scheme is non-dissipative and it is stable if \v\ < 1. (see the additional 
discussions given in Sec. 3). 

(f) For the dissipative extensions [5] of the non-dissipative forward marching form of the a scheme, only the 
less stringent local flux conservation condition Eq. (2.18) is imposed. Because CE(j, n), by definition, 
is a CCE, it follows that the solution to such an extension satisfies the following less stringent global 
flux conservation relation: the total flux of h* leaving the boundary of any space-time region that is the 
union of any combination of CCEs vanishes. Moreover, because Eq. (2.18) is equivalent to Eq. (2.12), 
for each of these extensions, it” is still evaluated using Eq. (2.12) while (%)" is evaluated using an 
equation different from Eq. (2.13). 

Next we will show that, in a sense to be defined shortly, any solution of the a scheme automatically 
satisfies the classical leapfrog scheme. Consider any (j. n + 1) € Hi. Then because (j,n+ 1) € fii -O- 
(j + 1, n) £ fli (j — l,n) £ Hi •<=>■ (j, ?i+ 1) £ fi i, Eqs. (2.12), (2.13), (2.16), and (2.17) imply that 


= 


»r' = 


{(1 - v) [u 

- (1 + V)u s ]j +1 + (1 + U) [U+ (1 - '')Wr]''.i}- 

(j) n + 1) € 

(2.19) 

.)«+! = 1/ 
■a 2 l 

[«-(! + v)ux\ n J+ 1 - [u + (1 - 

(j) n + 1) G Hi 

(2.20) 

{(1 + A) [u 

- (1 - v)u s ]j +1 + (1 - v) [u + (1 + 

(j, n - 1) £ Hi 

(2.21) 

,) n ~ 1 = -{ 
A 2 1 

[u - (1 - v)Ux] n j+ 1 - [u + (1 + z ')«*:]"_! j, 

(j,n- 1) £ Hi 

(2.22) 


and 


respectively. By subtracting Eqs. (2.21) and (2.22) from Eqs. (2.19) and (2.20), respectively, and using the 
definition v = aAt/Ax , one has two leapfrog schemes in which the mesh variables u" and (u^)” are decoupled, 
i.e., 

-A l — + a ^±L i_L = o, 


and 


2a t 

(u s )] +1 - (u,)] 
2a t 


2ax 
2ax 


= 0 , 


( ji n + 1) £ Hi 

( j , n + 1) £ Hi 


(2.23) 

(2.24) 


Hereafter, the system of equations defined by Eqs. (2.23) and (2.24) over all ( j,n + 1) £ fli will be referred 
to as the leapfrog form of the a scheme. Note that: (i) the a scheme implies the leapfrog form but not vice 
versa; and (ii) given the initial data at two consecutive time levels, the leapfrog form can be used for both 
forward and backward time marching. 

Let Uj and (u $ )°, j = ±1,±3, . . ., be given. Then Eqs. (2.12) and (2.13) can be used to define uniquely 

a set of it" and (%)" where ( j , n) £ Hi and n = 1.2. 3; On the other hand, with the aid of (i) the given 

initial data Uj and (u s )j, j = ±1, ±3, . . ., and (ii) the values of 1 1 ] and (w $ )j, j = 0, ±2, ±4, . . ., determined 
above, the leapfrog schemes Eqs. (2.23) and (2.24) can be used to define uniquely a new set of it” and (u s )" 
where (j,n) £ and n = 2,3,4, .. .. According to the result established in the last paragraph, for each 
(j, n ) £ Hi where n = 2, 3, 4, . . ., the values of w” and (u s .)" in the new set, respectively, should be identical 
to those in the set defined using the forward marching form of the a scheme. This conclusion has been 
verified in numerical experiments using periodic initial data (refer to the code listed in the appendix). 

This section is concluded with the following remarks: 

(a) Because, at each ( j , n) £ Hi, the a scheme is formed by two coupled equations involving two independent 
mesh variables it” and (%)", it is consistent with a pair of partial differential equations (PDEs) with 
one of them being Eq. (1.1) [1,11]. 

(b) The leapfrog scheme is a two-way marching scheme, i.e., given the data at any two consecutive time 
levels, Eq. (2.23) can be used for both forward and backward time marching. 
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(c) The leapfrog scheme Eq. (2.23) can be derived using a finite-difference [p.100, 73] or a finite- volume 
approach. In fact, for any (j, n) G fii, Eq. (2.23) can be cast into the local conservation form 

< f h*-ds = 0, (j + l,n)efii (2.25) 

JS(FV(j+ l,n)) 

where (i) F V(j + l,n) is the region ABCD depicted in Fig. 3; and (ii) the average flux vectors h* at 
AB, BC, CD, and DA are taking to be (au" +1 , u” +1 ), (. au'j_ 1 ,u'j _ 1 ), (au™ _1 , u™ -1 ), and (au” +1 , u” +1 ), 
respectively. Thus, for the leapfrog scheme, hereafter FV (j + l,n) is referred to as the conservation 
element (CE) assigned to the mesh point ( j + l,n). 

(d) Let (i) 

fin d = G fii and n = ±1, ±3, . . .} (2.26) 

(ii) 

$ii 2 d = {(j,n)\(j,n) G fii and n = 0,±2,±4, . . .} (2.27) 

def 

and (iii) En (Ei 2 ) = the set of the CEs of the leapfrog scheme assigned to all (j, n) G fin (fii 2 ). 
Then (i) fi i = fin U fii 2 , and (ii) the CEs in En and Ei 2 , separately , are non-overlapping and fill the 
entire space-time E 2 . Thus, by using Eq. (2.25) and the same argument used to establish the global flux 
conservation relation for the a scheme, one can establish two separate global flux conservation relations 
for the leapfrog scheme, i.e. , the total flux of h* leaving the boundary of any space-time region that is 
the union of any combination of CEs G On (Oi 2 ) vanishes. 

(e) The leapfrog scheme shares with the a scheme three key nontraditional features, i.e., they both are 
explicit two-way marching schemes, have space-time staggered stencils and their interface fluxes are 
evaluated without using any interpolation or extrapolation technique. However, due to its less sophisti- 
cated structure, it is much more difficult (if it is possible at all) to construct space-time flux-conserving 
dissipative extensions of the leapfrog scheme than those of the a scheme. Because (i) stable nonlinear 
simulations can only be achieved using dissipative extensions and (ii) a two level scheme is easier to 
handle than a three-level scheme, CESE development (which is based on the a scheme) is much simpler 
and more robust than that based on the leapfrog scheme. 

(f) By using Eq. (2.9), one can show that Eq. (2.24) is still valid if each (u$)", (j, n) G fii, is replaced by its 
mesh-size independent version (u x ) r t. Because the leapfrog scheme is second-order in accuracy for both 
space and time, one can infer from Eq. (2.23) and the mesh-size independent version of Eq. (2.24) that 
the a scheme is second-order in accuracy for both it” and (u^)™. In fact, this conclusion was proved 
rigorously as Eqs. (6.32) and (6.33) in [1], and can also be verified numerically using the code listed in 
Appendix. 

(g) In the above, the a and leapfrog schemes are defined using only the mesh points G O i. Independently, 
they can also be defined using only the mesh points G fi 2 where 

f h = {( j , n)|j, n = 0, ±1, ±2, ±3, . . . , and (j + n ) is an even integer} (2.28) 

For the current ID case where a structured mesh is used, both the a scheme and the leapfrog scheme 
defined over Hi are completely decoupled from those defined over fi 2 . Thus, there is no practical reason 
to carry out computations using the two decoupled schemes simultaneously . Nevertheless, the leapfrog 
scheme is traditionally referred but unjustifiably to that formed by the two decoupled leapfrog schemes. 
Moreover, through a twisted logic, the leapfrog scheme is falsely blamed for its solution becoming 
decoupled after a long time marching. 

To simplify mathematical manipulation, hereafter we will exclusively deal with the “dual” a scheme 
and the “dual” leapfrog scheme which , respectively, are formed from two decoupled a schemes and two 
decoupled leapfrog schemes. As such these dual schemes are defined over 

fi = fliU0 2 = {(j,n)\j,n = 0,±1,±2,...} (2.29) 
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Moreover, from now on it is to be understood that (i) the a scheme and the leapfrog scheme refer to 
the dual schemes which are defined over Q ; and (ii) the “dual” version of a previously defined system of 
equations will be identified using the same equation number of the original version with an extra symbol 
“d” being attached at the end. 

3. von Neumann Analysis 

As a preliminary to later developments, the results of the von Neumann analysis for several non- 
dissipative (i.e., neutrally stable) solvers for Eq. (1.1) will be studied in this section. In particular, we 
will point out key properties shared by the amplification matrices (factors) of these schemes. 

The Wendroff scheme [p.503, 72], i.e., 

(1 + v){u™+l - u ™ ) + (1 - v){u™ +l - u" +1 ) = 0, (j, n) € Cl (3.1) 

is a two-level and two-way marching implicit non-dissipative solver of Eq. (1.1). Let 

u] = [g(iy,9)] n e^, (j,n) G Q (3.2) 


be a solution to Eq. (3.1) where (i) i = \f—l, (ii) 9 (— oo < 9 < oo) is a phase angle, and (iii) g{v,9) (^ 0) 
is a complex function of v and 9 (note: [ g{v , 9)} n is undefined for any n < 0 if g(y : 9) = 0). By substituting 
Eq. (3.2) into Eq. (3.1), one has 


[(! + v)(ge l9 - 1) + (1 - v)(g - e l6 )} g n = 0, n = 0, ±1, ±2, . . . , 

Because g° = 1, Eq. (3.3) 

[(1 + u)e ie + (1- V )\g = (i- v)e ie + (1 + v) 


(3.3) 


(3.4) 


Dividing the expressions on both sides of Eq. (3.4) by 


e i 0 / 2 anf ] us j n g b as i c trigonometry, one has 


def cos(0/2) — ivs\n{9/2) 
cos(0/2) + ivsvn(6/2) 


(cos(0/2) ± ivsm{0/2) ^ 0) (3.5) 


Note that: (i) cos(0/2) + ivsm(Q/2) = 0 cos(0/2) — ivs\n{9/2) = 0 v = 0 and 9 = ±7r, ±37 t, ±57t, . . ., 

i.e., the domain of g w {v,9) is 


D{g w ) = {(u, 9)\v ^ 0 with any real 9 , or v = 0 with any real 9 ^ ±i r, ±37r, ±57 t, } (3.5a) 

(ii) by definition, g w (is,9) is the amplification factor of the Wendroff scheme; and (iii) over D(g w ), we have 


g(v, ~9) ■ g{v, 9) = 1 (3.6) 

and 

g(v,9) = g{v,~9) (3.7) 

if g(v,9) = g w (v,9). Here, for any complex number c, its complex conjugate is denoted by c. Combining 
Eqs. (3.6) and (3.7), one concludes that, for all real 9, we have 


\g(v,9)\ = l (3.8) 

if g(y,9) = g w {v,9) and v ^ 0. Thus, practically, the Wendroff scheme is unconditionally neutrally stable. 

At this juncture, note that the property Eq. (3.7) is a result of the fact that the Wendroff scheme is 
a linear solver of Eq. (1.1) which has real constant coefficients and a single mesh variable per mesh point. 
This can be seen by the following arguments: By taking the complex conjugate of Eq. (3.1), it is seen that 
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the complex conjugate of any solution to the solver must also be a solution. Thus, if Eq. (3.2) is a solution, 
so is 


g(v, 0 ) 


?~ijO 


(. j , n) G n 


(3.9) 


Because (i) for any (v, 6) G D{g w ), g{v,9) is uniquely determined by the requirement that Eq. (3.2) is a 
solution to Eq. (3.1), and (ii) (v,9) G D{g w ) <G> (v,—9) G D{g w ), a comparison between Eq. (3.2) and the 
equation resulting from replacing 9 with —9 in Eq. (3.9) reveals that, for any (i ;,9) G D{g w ), 


g{v,-9) = g(v,9) 


(3.7 a) 


Again because (v,9) G D{g w ) <G> {v,—9) G D{g w ), Eqs. (3.7) and (3.7a) are equivalent. QED. As such, the 
property Eq. (3.7) is shared by any solver of Eq. (1.1) which shares the same characteristics of the Wendroff 
scheme defined above. Because this solver can even be a dissipative solvers such as the Lax scheme [p.97, 73] 
the key factor contributing to the neutral stability of the Wendroff scheme is Eq. (3.6) rather than Eq. (3.7). 

Another two- level and two-way marching implicit non-dissipative solver of Eq. (1.1) is the Crank- 
Nicolson scheme [p.504, 72], i.e., 

u] +1 - U 1 + + «"+i - <-i - «"-i) = 0, (j, n) G Q (3.10) 

Through a process similar to that presented above, it can be shown that the amplification factor for the 
Crank-Nicolson scheme is 


9cn{v 1 0) 


def 1 — {iv/ 2) sin# 
1 + {iv/ 2) sin#’ 


—oo < v,9 < +oo 


(3.11) 


It can be shown easily that, for any real v and 0, Eqs. (3.6)-(3.8) are satisfied if g{v,9) = g C n{v,9). Thus 
the Crank-Nicolson scheme is unconditionally neutrally stable. 

There is no classical two-level neutrally stable explicit solver of Eq. (1.1). In fact, among the classical 
schemes, the three-level leapfrog scheme is the only well-known neutrally stable explicit scheme. To study 
its amplification factors, first we will cast it into an equivalent two-level forward marching matrix form. As 
a preliminary, let 

P d ='(; J), and «?(; (3.12) 

and 

<f(j, n) d = ’ (i> n ) G O (3.13) 

where (/>" and </>", {j,n) G O are arbitrary real constants. Then (i) 

P = P~ 1 (3.14) 


where P 1 is the inverse of P; and (ii) the system of matrix equations 


n+l) = P ${j, n) + H{v) ${j - 1, n) - (j>{j + 1, n) 


{j,n)efl (3.15) 


is equivalent to 
and 


0" +1 = $ + v [ct>U - <% +1 ] , (j, n) G O 

{j, n) G O 


P] +1 = 4%, 


(3.16) 

(3.17) 
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From the above observation (ii), one concludes that Eqs. (3.13) and (3.15) are equivalent to 


n) = ( p- 1 ) » (j, n) € (3.18) 

and 

0” +1 = <T 1 + " l >?- 1 - <A "+1 ] , (i, n) e o (3.19) 

By comparing Eq. (3.19) with Eq. (2.23d) and using Eq. (3.18), it is seen that Eq. (2.23d) is equivalent to 
Eq. (3.15) if </>" is taken to be u ”, i.e., if 


${j, n) d = ^ i ^ , (j, n) € O (3.20) 

Hereafter, assuming Eq. (3.20), Eq. (3.15) is referred to as the forward marching form of the leapfrog scheme. 
For any real v and 9, let A(v 1 9) be a 2 x 2 nonsingular complex matrix function of v and 9 such that 

${j, n) = e ljB [A(u, 9)] n a, ( j , n) G O (3.21) 

is a solution to Eq. (3.15) for all possible complex constant 2x1 column matrices a (note: \A(v, 9)] n is 
undefined for any n < 0 if A(u, 9) is singular). By substituting Eq. (3.21) into Eq. (3.15), one has 

[A(v,9)-P-H(v)(e~ iB -e ie )\[A(v,9)} n a = 0, n = 0, ±1, ±2, . . . (3.22) 

Because (i) [A(i/, 0)]° = I where I is the 2x2 identity matrix, and (ii) a can be any complex constant 2x1 
column matrix, with the aid of Eq. (3.12), one concludes that Eq. (3.22) <t=> 

A(v,6) = P + H(v)(e~ ie - e i0 ) = ^~ 2w ^ md ^oo<^,0<+oo (3.23) 

According to Eq. (3.23), det [A(i/, 0)] = — 1. Thus A(u. 9) is indeed nonsingular over its domain. By definition, 
A(v. 9) is referred to as the amplification matrix of the leapfrog scheme with respect to the form Eq. (3.15). 
It can be shown easily that 


A(v, -0)PA(v, 9)P~ X = /, -oo < v, 9 < +oo (3.24) 

and 

A(v, 9) = A(u, — 9), —oo<i , 1 9<+oo (3.25) 

Here, for any matrix M, its complex conjugate is denoted by M. Combining Eqs. (3.24) and (3.25), one has 

A(v, 9) PA(v, 9)P^ 1 = /, — oo < i/, 0 < +oo (3.26) 

Eqs. (3.24) — (3.26) are the current counterparts to Eqs. (3.6)-(3.8), respectively. 

By definition, the two eigenvalues of A(v. 9) are referred to as the amplification factors of the leapfrog 
scheme Eq. (2.23d). Let these factors be denoted by \+(v, 9) and A_ (v, 9), respectively. Then the eigenvalues 
of A(v,9) are A +(v,9) and \-(v, 9), respectively. On the other hand, because PA(v, 9)P~ 1 is similar to 
A(v, 9) (i.e., these two matrices are related by a similarity transformation), the eigenvalues of PA(y : 9)P ~ 1 
are also \+{v, 9) and A _(i/, 0), respectively [76]. Moreover, we have 


det 


A(u, 9) ■ det [PA(n, 9)P 1 = det (I) = 1, 


—oo < is, 9 < +oo 


(3.27) 
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which is a result of Eq. (3.26) and the matrix theorem: det (AB) = det(A) • det(H) for any N x N matrices 
A and B. By using the above observations and another matrix theorem, i.e., the determinant of a square 
matrix is equal to the product of the eigenvalues of the same matrix, one now arrives at the conclusion that 

A +(is, 9) ■ \-{v, 9) ■ A +(v, 6) ■ A _ (v, 9) = 1, — oo < u,9 < +oo 


i.e., 


| A + 0)| • | A_ (zv, 0)| = 1, — oo < v, 9 < +oo (3.28) 

For any given v, stability of the leapfrog scheme requires that 

|A + (z',0)| < 1 and |A_(;/,0)| < 1, — oo <9 < +oo (3.29) 

Thus Eq. (3.28) implies that, for any given v, the leapfrog scheme must be neutrally stable, i.e., 

\A+(v, 9)\ = |A_(^, 0)| = 1, —oo <9 < +oo (3.30) 

if it is stable. As such, Eq. (3.26) does not imply neutral stability of the leapfrog scheme. However, it does 
imply that the scheme can only be neutrally stable (i.e., non-dissipative) if it is stable. Here we have reached 
this conclusion without using the explicit form of A+(v, 9) and A-(v, 9), i.e., 


A±(u,9) = — iu sin 9 E \/l — v 2 sin 2 9 


(3.31) 


In fact, by using Eqs. (3.23) and (3.31) along with a rigorous definition of stability given in [64], it can be 
shown that the leapfrog scheme is (i) neutrally stable if v 2 < 1; (ii) linearly unstable when v 2 = 1; and (iii) 
exponentially unstable if v 2 > 1. Note that the linear instability at v 2 = 1 is due to the fact that A(v, 9) 
becomes defective when v 2 = sin 2 9=1 (see pp. 6-9 in [64]). 

Note that, except for some simple ID cases, explicit forms of amplifications factors generally are un- 
available. However, as explain earlier, the matrix relation such as Eq. (3.26) may be used to draw important 
conclusions about possible neutral stability of a scheme even if no explicit amplification factors are available. 

Next we will show that: (i) the leapfrog scheme Eq. (2.23d) can be cast into different equivalent two- 
level matrix forms; and (ii) the amplification matrices associated with different matrix forms are related by 
similarity transformations and therefore they share the same eigenvalues, i.e., the same amplification factors. 

Let 




x def 

n) = 


u i 

i /" -1 

u j± i 


(j, n) G fl 


and 


T def 
-/ + = 


tt / \ def / V 1 

H+ M =00 


def 


1 

0 A 

def ( 0 

0 A 

o 

O’ 

i— « 

1 

II 

O 

l 


tt ( \ def / — V 0 

H -M =10 


R(9) = I + + I_ e ie , -oo < 9 < +oo 

Eqs. (3.33) and (3.35) imply that 

( I +) 2 = I +, (/-) 2 = /_, /_/+ = /+/_= 0, /_+/+=/ 

and 


(3.32) 

(3.33) 

(3.34) 

(3.35) 

(3.36) 


(. R(9))~ 1 = R(—9), -oo < 9 < +oo (3.37) 

where 0 is the 2x2 null matrix. Moreover, by using Eqs. (3.20) and (3.32)-(3.34), one can shown that: (i) 


<f>±(j, n) = I+(j>{j, n) + I-$(j ± 1, n), 


(j, n ) G Q, 


(3.38) 
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and (ii) Eq. (2.23d) is equivalent to each of the following two forms: 


$+(j,n + 1) = H-(i/)$ + (j + l,n) + H+(i/)$+(j - l,n), (j,n)efl (3.39) 

and 

1) = [H^v)]* + \,n) + [H+(v)] t -\,n), (. j,n ) G O (3.40) 

where [ H±(u)] t is the transpose of H±{v). 

Next, by substituting Eq. (3.21) into Eq. (3.38) and using Eq. (3.35), one concludes that Eq. (3.21) can 
be cast into the following equivalent forms: 

$±{j, n) = e l36 R{±0) [A(v, 9)} n a, (j, n) G (3.41) 

According to Eqs. (3.37), R(=f6)R(±9) = I. Thus Eq. (3.41) can be recast as 

$±{j,n) = e lj6 [A±(v,9)] n a±, (j,n)eQ (3.42) 

where 

A±{v,9) = f R(±9)A(u,9)R(t 9) = R{±9)A(v,9) [i?(±6»)] _1 (3.43) 

and 

a± d = R(±9)a (3.44) 

Note that: (i) Eq. (3.21) is a solution to Eq. (3.15) for all possible complex constant 2x1 column matrices a 
<t=> Eq. (3.23); (ii) Eqs. (3.15), (3.39), and (3.40) are equivalent forms of the same leapfrog scheme Eq. (2.23d), 
and (iii) non-singularity of R(±6) implies that a± can be any complex constant 2x1 column matrix just 
like a. As such, a comparison between Eqs. (3.21) and (3.42) reveals that A±(v,9') defined by Eqs. (3.23), 
(3.35), and (3.43), i.e., 


A+{v, 9) 


— 2ivsm 9 

A6 



and A _ ( v , 9) 


— 2ivsin9 e ie \ 
e~ i9 0 ) ’ 


—oo < is, 9 < +oo 


(3.45) 


must be the unique solutions to Eqs. (3.39) and (3.40), respectively. This fact can also be verified directly 
by substituting Eqs. (3.42) into Eqs. (3.39) and (3.40), respectively. Thus A + {v,9) and A_(is, 9) defined 
by Eq. (3.45) are the amplification matrices of the leapfrog scheme with respect to the forms Eqs. (3.39) 
and (3.40), respectively. Moreover, (i) because A(u, 9), A + (v,9), and A_(v, 9) are related by similarity 
transformations Eq. (3.43), they have the same eigenvalues (i.e., amplification factors) A + (z/, 9) and A_(z^, 9)\ 
and (ii) by using Eqs. (3.45), one can show that Eqs. (3.24)-(3.26) are still valid if A(u, ±9) is replaced by 
A + {y,±9) or A_(v,±9). 

Next we consider the a scheme. Let 


-*( . \ def 

q[j,n) = 


U J 

(%)j 


and 


, x def 1 / 1 + V 1 - V 2 

e+ ( o = 2 ( — i -(i-o 


(. j , n) G ft 


Q_(„) d £ 


(3.46) 


(3.47) 


1 -(1 + ^) 

Then the forward marching form Eqs. (2.12d) and (2.13d) of the a scheme can be cast into the matrix form: 


q(j,n) = Q-{v)q(j + l,n- 1) + Q+(v)q(j - l,n- 1), (j,n) G ft 

Let G(z/, 9) be a 2 x 2 nonsingular complex matrix function of v and 9 such that 

q{j , n) = e lje [G(v, 9)] n b, ( j , n) G O 


(3.48) 


(3.49) 
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is a solution to Eq. (3.48) for all possible complex constant 2x1 column matrices b. By substituting Eq. (3.49) 
into Eq. (3.48), one has 


\G(v, 9) — e ie Q-(y) — e~ ie Q + (v)] [G(v, 9)] n b= 0, n = 0, ±1, ±2, . . . (3.50) 


Because (i) \G(y, 0)]° = I, and (ii) b can be any complex constant 2x1 column matrix, with the aid of 
Eq. (3.47), one concludes that Eq. (3.50) <t=> 


G[v,&) = e ie Q-{v) + e~ w Q + {v) 


f cos 9 — iv sin 9 
y i sin 9 


—i(l — v 2 )sm0 \ 
— (cos 6 + ivsui9) ) ’ 


—oo < v, 9 < +oo (3.51) 


According to Eq. (3.51), det[G(+ 0)] = —1. Thus G(y, 9) is indeed nonsingular over its domain. By definition, 
G(v, 9) is referred to as the amplification matrix of the forward marching form of the a scheme. It can be 
shown that the eigenvalues of G(v,0) are A +(v,0) and A ~{v,9) defined in Eq. (3.31), i.e., the amplification 
factors of the forward marching form of the a scheme are identical to those of the leapfrog scheme. 

Next let 


jj def / 1 0 

L = 1,0 -1 

) 

(3.52) 

and observe that 

U = U~ l 


(3.53) 

Thus Eq. (3.51) implies that 



G(v,-9)UG{v,6)U~ l = 1, 

—oo <v,9< +oo 

(3.54) 

and 

G(v,e)=G(v,-9), 

Combining Eqs. (3.54) and (3.55), one has 

-oo <v,9< +oo 

(3.55) 

G(v,9) UG(v, 9)U~ 1 = I, 

—oo <v,9 < +oo 

(3.56) 


Eqs. (3.54)-(3.56) are the current counterparts to Eqs. (3.24)-(3.26), respectively. 

Eqs. (3.26) and (3.56) are identical in their forms. Thus, without using the explicit forms of A + (y,9) 
and A_(i/, 9), the same argument that was given following Eq. (3.26) can be invoked to show that, for any 
given v, the forward marching form of the a scheme must be neutrally stable if it is stable. In fact, by using 
Eqs. (3.51) and (3.31) along with a rigorous definition of stability given in [64], it can be shown that the 
forward marching form of the a scheme is (i) neutrally stable if v 2 < 1; (ii) linearly unstable when v 2 = 1; 
and (iii) exponentially unstable if v 2 > 1. Note that the linear instability at v 2 = 1 is due to the fact that 
G(y, 9) becomes defective when v 2 = sin 2 9 = 1. 

At this juncture, note that Eqs. (3.25) and (3.55), like Eq. (3.7), follow directly from a rather common 
property that Eqs. (3.15) and (3.48) are linear homogeneous systems of equations with real constant coeffi- 
cients. Thus, one concludes that Eqs. (3.6), (3.24), and (3.54) are the deciding contributing factors to the 
non-dissipatedness of the Wendroff, Crank-Nicolson, leapfrog and a schemes. In Sec. 4, we will show that 
these latter equations are the natural results of a unique space-time invariance property shared by the above 
non-dissipative schemes. To pave the way, this section is ended with the following remarks: 

(a) The leapfrog scheme Eq. (2.23d) can also be cast into a two-level backward marching matrix form, i.e., 


V’lb n ) = n + 1) + H (v) ip(j + 1, n + 1) - ip(j - 1, n + 1) 


(j, n) G O 


(3.57) 


where 


Wj,n) = 




(j, n) G 


(3.58) 
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(b) The backward marching form Eqs. (2.16d) and (2.17d) of the a scheme can be cast into the matrix form: 


where 


and 


n) 


Q+(v)q(j + l,n+ 1) + Q-(v)q(j - l,n+ 1), (j,n) G O 


(3.59) 

(3.60) 

(3.61) 


4. Space-time inversion invariance and non-dissipatedness 

Note that Eq. (1.1) is equivalent to 


where 


du du 

df +a d^ = 


x' = f 2x c — x and t! = f 2 t c — t 


(4.1) 

(4.2) 


with x c and t c being any pair of given real constants. Thus Eq. (1.1) is mapped into the same equation 
under the mapping 

(. x,t ) <-> {x' ,t') (4-3) 


Because Eq. (4.2) <G> x + x' = 2x c and t + t' = 2 t c , the transformation defined in Eq. (4.2) represents a 
space-time inversion (i.e., a spatial reflection followed by a time reversal) with respect to the point ( x c ,t c ). 
As such Eq. (1.1) is said to possess the space-time inversion (STI) invariance. 

As a result of the above invariance property, solution to Eq. (1.1) also possesses a similar invariant 
property. Let 

u=u 0 (x,t) (4-4) 


be a solution to Eq. (1.1). Then because (i) u = u 0 (x,t) is a solution to Eq. (1.1) <G> u = u 0 (x',t') is 
a solution to Eq. (4.1), and (ii) Eqs. (1.1) and (4.1) are equivalent if x ’ and t' are related to x and t by 
Eq. (4.2), one concludes that 

u = u 0 {2x c — x, 2 t c — t) (4-5) 


must also be a solution of Eq. (1.1). In fact, the general solution to Eq. (1.1) is given by 


u = F(x — at) 


(4.6) 


where F(s) is an arbitrary differentiable function of s. It is easy to see that u = F({ 2x c — x) — a(2t c — t)) is 
also a solution to Eq. (1.1) 

Numerical schemes generally are constructed without considering the invariant properties of the physical 
equations they model. Thus numerical solutions generally do not possess the same invariant properties of 
physical solutions. As an example, the general solution Eq. (4.6) represents non-dissipative data propagation 
along the characteristic lines defined by dx/dt = a. However, numerical solutions generated by dissipative 
analogies of Eq. (1.1) are dissipative. Using the numerical analogies of Eq. (1.1) defined earlier as examples, 
it will be shown in this section that non-dissipatedness of a scheme is deeply related to its STI invariance 
property. 

Let j c and n c denote any pair of given fixed integers. Then the mapping 


(j, n) <-> (2 j c - j, 2 n c - n), (j, n) G ft (4.7) 

represents a space-time inversion mapping of mesh points (j, n) G U with respect to ( j c , n c ). As an example, 
under the mapping Eq. (4.7) with j c — n c = 0, points A , H, C, and D depicted in Fig. 4 map onto points 
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A', B ', C' , and D’ , respectively, and vice versa. Thus a system of equations, such as Eq. (3.1), is said to be 
STI invariant if the system maps onto itself under the one-to-one mapping 

^ (j,n)£fl (4.8) 

In other words, a STI invariant system is a system which will transform back to the original system if, for 
each (j,n) £ Cl, Uj is replaced by u^Zj 1 where j c and n c are any pair of given fixed integers. 

Let j' j — j c and n! = f n — n c . Then (i) j = j c and n = n c <t=> j' = n' = 0 and (ii) (j, n) G Cl 
<G> ( j',n ') £ Cl. As a result, without any loss of generality, in the following developments we can and will 
consider only the special case of Eq. (4.8) with j c = n c = 0, i.e., 


( j,n)£Cl (4.9) 

First we shall establish the STI invariance of the Wenclroff scheme. Under the mapping Eq. (4.9), 
Eq. (3.1) is mapped onto 

(1 + v) - uzfj + (1 - v) - «I" j+1) ) = 0, ( j , n) G Cl (4.10) 

(note: Eq. (4.10) is not the equation that we would obtain from Eq. (3.1) by simply replacing the symbols 
“ j ” and “n” in it with —j and — n, respectively). Let 


j* = -(j + l) and n* = — (n + 1), (j,n) G Cl 


(4.11) 


Then, with the aid of Eq. (4.11) and the fact that (j* ,n*) G 0 <t=> ( j,n ) G Cl, Eq. (4.10) can be rearranged as 

(1 + v) (ufXl - uf) + (1 - v) (uf, +1 - uf, + 1 ) = 0, C r,n *) G Cl (4.12) 

By comparing Eq. (4.12) with Eq. (3.1), one can see that the system of equations defined by Eq. (4.12) (and 
therefore that defined by Eq. (4.10)) is identical to that defined by Eq. (3.1). QED. 

To gain a deeper insight about STI invariance, consider the two component equations of Eq. (3.1) with 
{j,n) = (1,1) and (J,n) = (-2,-2). They are 

(1 + v){u% - «J) + (1 - Z ')(«? - u\) = 0 (4.13) 

and 

(1 + v)(u_ 2 — u_i) + (1 — v)(u_\ — u_ 2 ) = 0 (4-14) 

respectively. Obviously Eq. (4.14) is the image of Eq. (4.13) under the mapping Eq. (4.9) and vice versa. 
This represents a concrete example of the fact that the image of a component equation of Eq. (3.1) under 
the mapping Eq. (4.9) is also a component equation of Eq. (3.1), and vice versa. 

Moreover, (i) the mesh points involved in Eq. (4.13) are points A, B , C, and D depicted in Fig. 4 
while those involved in Eq. (4.14) are points A' , B' , C' , and D'\ and (ii) under the mapping Eq. (4.7) with 
jc = n c = 0, i.e., 

(j,n) <-> {-j,~n), (. j,n)£Cl (4.15) 

points A' , B ' , C' , and D' are the images of A, B, C, and D, respectively, and vice versa. Because, for a 
linear constant-coefficient solver of Eq. (1.1) with one mesh variable and one equation per mesh point, the 
stencil of all its component equations have the same configuration, the above observation implies that such 
a scheme cannot be STI invariant unless it meets the requirement that its stencil configuration is unchanged 
under the space-time inversion mapping Eq. (4.9). Obviously this requirement is met by the non-dissipative 
Wendroff, Crank-Nicolson, and leapfrog schemes. 
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On the other hand, any two-level explicit linear constant-coefficient (forward marching) solver ofEq. (1.1) 
with one mesh variable and one equation per mesh point can never meet the above requirement. By its ex- 
plicit nature, the stencil of such a scheme can contain only one point at the upper time level. Whereas 
stability and other considerations demand that the stencil contains at least two points at the lower time 
level. As such, the configuration of its stencil will change under the mapping Eq. (4.15). In turn, this implies 
that such a scheme cannot be STI invariant. In view of the strong tie to be established between STI invariance 
and non-dissipatedness, this explains why none of such two-level explicit schemes is non- dissipative. 

Next we will establish a relation between STI invariance and non-dissipatedness. Let s(j,n), ( j,n .) G fi, 
be a function of (j, n ) such that 

u] = s(j, n), (j, n) G fi (4.16) 

is a solution to Eq. (3.1). Then, by substituting Eq. (4.16) into Eq. (3.1), one obtains the system of identities: 

(! + v) ( s(j + 1, n + 1) - s(j, n)) + (1 - v) ( s(j , n + 1) - s(j + 1, n)) = 0, ( j , n) G ft (4.17) 

On the other hand, because Eq. (4.10) <G> Eq. (3.1), by substituting Eq. (4.16) into Eq. (4.10), one concludes 
that the system of identities 

(1 + v) (s(-(j + l),-(n+ 1)) - s(-j,-n)) + (l-u) (s(-j,-(n+ 1)) - s(-(j + 1 ),-n)) = 0, (j,n) G ft 

(4.18) 

is identical to that defined in Eq. (4.17). As a concrete example, the component identities in Eq. (4.17) with 
(j,n) = (1,1) and (j,n) = (-2,-2) are 

(1 + v) (s(2, 2) - a( 1, 1)) + (1 - v) (s(l, 2) - s( 2, 1)) = 0 (4.19) 

and 

(1 + v) (s(— 2, -2) - s(- 1, -1)) + (1 - v) (s(- 1, -2) - s(— 2, -1)) = 0 (4.20) 

respectively. Whereas the component identities in Eq. (4.18) with (j, n) = (1,1) and (j, n) = (—2,-2) are 

Eqs. (4.20) and (4.19), respectively. 

Moreover, a direct comparison of Eqs. (3.1) and (4.18) reveals that 

Uj = s(-j , -n), ( j , n) G ft (4.21) 

represents another solution to the system of equations Eq. (3.1), i.e., Eq. (4.16) is a solution to Eq. (3.1) <G> 
Eq. (4.21) is a solution to Eq. (3.1). As an example, by substituting Eq. (4.16) into Eq. (4.13) and (4.14), 

one has Eqs. (4.19) and (4.20), respectively. Whereas, by substituting Eq. (4.21) into Eqs. (4.13) and (4.14), 

one has Eqs. (4.20) and (4.19), respectively. In the following, the key property Eq. (3.6) of the Wendroff 
scheme will be derived using the fact just established above without using the explicit form of g w (v, 9) given 
in Eq. (3.5). 

Let Eq. (3.2) be a solution to Eq. (3.1) for all (v,9) G D(g w ). Then g(is,9) must take the unique form 
defined in Eq. (3.5). Moreover, the fact established in the last paragraph implies that, for all (v, 9) G D(g w ), 

u] = [g(u,9)]- n e- ije = [(g(v, S))” 1 ]" e~ ij \ (j ,n) G (4.22) 

must also be a solution to Eq. (3.1). Because (v,Q) G D(g w ) <G> (u, —9) G D(g w ), by replacing 9 in Eq. (4.22) 
with —0, one concludes that, for all (n. 9) G D(g w ), 

u] = l(g(v-9))- l ] n ji 9 , ( j,n ) G O (4.23) 

must also be a solution of Eq. (3.1). By comparing Eq. (3.2) and (4.23) and recalling the uniqueness of 
g(v,9), one concludes that g(v,9) = (g(v 1 —9))~ 1 , i.e., Eq. (3.6). QED. 
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Using arguments similar to those just presented, one can easily show that the Crank-Nicolson scheme is 
also STI invariant and possesses all other properties established above for the Wenclroff scheme. In particular, 
it can be shown that the fact that g(y,9) = g C n(y,9) satisfies Eq. (3.6) is a result of its STI invariance. 

Similarly, one can easily establish STI invariance for the leapfrog scheme Eq. (2.23d) and show that 
Eq. (4.16) is a solution to Eq. (2.23d) <+> Eq. (4.21) is a solution to the same equation. However, as will be 
shown later, it is trickier to establish the key relation Eq. (3.24). 

As a preliminary, first we offer the following remarks: 

(a) Let Eq. (4.16) be a solution to Eq. (2.23d). Then Eq. (3.20) implies that 

*’" ) = UIn-l))' < 4 ‘ 24 > 


is a solution to Eq. (3.15), an equivalent form of Eq. (2.23d). Moreover, because Eq. (4.16) is a solution 
to Eq. (2.23d) <+> Eq. (4.21) is a solution to Eq. (2.23d), one concludes that 




( s(-j, -n) \ = / s(-j, 1-n) \ 

\ s (—j, 1 — n) J 


(. j , n) G 


(4.25) 


must also be a solution to Eq. (3.15). 

(b) Let j 0 and n 0 be arbitrary given fixed integers. Then (j + j 0 ,n + n 0 ) eflo (j, n) G fi. Thus, by direct 
substitution, one can show that, for any real v and 9, Eq. (3.21) is a solution to Eq. (3.15) <+> 

${j, n) = e %>+M e [A{y, 9)} n+n ° a, ( j , n) G U (4.26) 


is a solution to the same equation. 

For any real v and (9, let Eq. (3.21) be a solution to Eq. (3.15). Then, as was shown earlier, A{y , 9) must 
be in the unique form defined in Eq. (3.23). Also we can assume that: 


( S (i,n-1)) = e y V(M)] n a, l?» G Q (4.27) 

By replacing the indices j and n in Eq. (4.27) with — j and 1 — n, respectively and using the fact that 
(j, n) G (— j , 1 — n) G ft, one concludes that Eq. (4.27) <t=> 

(s(-T(i -n)~ 1)) = ^[AW)] 1 -^, (j,n) G O (4.28) 

By using Eq. (4.28) and the above remark (a), one concludes that, for any real v and 9 , 

$(j, n) = Pe~ i30 [A{v, 9)] l ~ n a = e~ iie P[A(v, 9)} 1 ~ n a, ( j , n) G Q (4.29) 


must also be a solution to Eq. (3.15). 

To cast Eq. (4.29) into another form, let 


A(v, 9) = f PA{v, 9)P x , — oo < v, 9 < +oo 


(4.30) 


Then, with the aid of the relation 


[PA(v, 9)P~ 1 ] 1 = P[A^,9)]~ l p-\ 


— oo < v,9 < +oo 


one can show that 


P[A(i 


A(is,0) P, n = 0, ±1, ±2, . . . ; — oo < v, 9 < +oo 


(4.31) 


(4.32) 
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By using Eq. (4.32), Eq. (4.29) can be expressed as 


= e lj6 \A{v, 0) 


I 1 —n 


Pa, 


(j, n) G 


(4.33) 


Then, by replacing 0 in Eq. (4.33) with — 0, it is seen that, for any real v and 0, 

I 1 —n 


Pa, 


(. 3 , n) G 0 


= e l ° e |A(i/,-0)J 

must be a solution to Eq. (3.15). In turn, by using the above remark (b), One concludes that 
= e zj6 \A(v , -0)1 Pa = \{A{v, -0)) _1 j Pa, ( j,n ) G O 


(4.34) 


(4.35) 


must be a solution to Eq. (3.15). Because P is nonsingular, Pa can be any complex constant 2x1 column 
matrix just like a. Thus, by (i) comparing Eq. (4.35) with Eq. (3.21), (ii) using Eq. (4.30), and (iii) recalling 
the uniqueness of A{y, 0), one has 

A(u, 0) = [PA(v, — 0)P _1 ] 1 = P[A{v, — 0)]~ 1 P~ 1 , — oo < v, 0 < +oo (4.36) 


which, by using Eq. (3.14), can be shown to be equivalent to Eq. (3.24). QED. 

Next, we shall define and establish STI invariance for the a scheme. Note that: (i)du/dx «-> —du/dx 
under the mapping x <-> — x ; and (ii) (ux)™ is the numerical analogy of the value of du/dx at the mesh 
point {j, n) . Because of the above observations and the definition Eq. (2.9), the a scheme is said to be STI 
invariant if and only if the associated system of equations maps onto itself under the one-to-one mapping 



By using Eqs. (3.46) and (3.52), Eq. (4.37) can be expressed as 

q{j, n) <-> £/ q{-j, -n), ( j , n) G Q 


(4.37) 


(4.38) 


Under the mapping Eq. (4.38), the forward marching matrix form Eq. (3.48) is mapped onto 

Uq(-j, -n) = Q -{v)U q{— j - 1, -n + 1) + Q+(v)Uq{-j + 1,-n+l), ( j,n ) G tt (4.39) 

Let Eq. (4.39) be multiplied from left by the nonsingular matrix U. Then, with the aid of Eqs. (3.53), (3.60) 
and (3.61), one concludes that Eq. (4.39) is equivalent to 

= Q+{y)q(—j + 1) ~n + 1) + Q~{v)q{—j - 1,-n + l), (j,n) G U (4.40) 

Moreover, because (i) (j, n) G O <+> (— j, — n) G U, and (ii) by replacing the indices j and n with —j and 
— n, respectively, Eq. (3.59) will turn into Eq. (4.40), one concludes that the system of equations defined by 
Eq. (4.40) is identical to that defined by Eq. (3.59), another equivalent form of the a scheme. Thus the a 
scheme is STI invariant. QED. 

Next, we will establish Eq. (3.54) from the STI invariant property of the a scheme. For any real v and 
0, let Eq. (3.49) be a solution to Eq. (3.48). Then the STI invariance of the a scheme and Eq. (4.38) imply 
that, for any real v and 0, 

q(j, n) = Ue-W [G(v, 9)}~ n b = e~^U [G(v, 0)]” n b, ( j , n ) G Q (4.41) 
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must also be a solution to Eq. (3.48). By using a line of arguments similar to that used to obtain Eq. (4.36) 
from Eqs. (3.21) and (4.29), Eq. (3.49) and (4.41) can be used to show that 

G(v, 0) = [UG(y, — 6)U ~ l ] _1 = U[G(u, -oo < v, 9 < +oo (4.42) 

By using Eq. (3.53), it can be shown that Eq. (4.42) is equivalent to Eq. (3.54). QED. 

This section is concluded with the following remarks (see Fig. 5): 

(a) Let (j,n) = (2,2). Then (j,n) is point A. According to Eqs. (2.10d), (2. lid), (2.14d), and (2.15d), the 
basic stencils of the a scheme which involve point A are formed by the following four pairs of mesh points: 
(i) A and B , (ii) A and C, (iii) A and D , and A and E. The forward marching relations Eqs. (2.12d) 
and (2.13d) are derived using the basic relations involving the mesh variables of pairs (iii) and (iv), 
while the backward marching relations Eqs. (2.16d) and (2.17d) are derived using those involving the 
mesh variables of pairs (i) and (ii). 

(b) The images of points A , B, C, and D under the mapping Eq. (4.15) are points A', B', C' , and D', 
respectively. Obviously, the configuration of the image of any stencil referred to in the above item (a) 
is identical to that of the original stencil. 

(c) From the observations given above, it becomes clear that a two-level explicit solver of Eq. (1.1) cannot 
be STI invariant unless it has at least two independent mesh variables assigned to each mesh point. 

5. Conclusions and Discussions 

A thorough and rigorous discussion of the core ideas and the core a scheme of the CESE method has 
been presented. These simple core ideas provide the CESE method a foundation which is solid in physics and 
yet simple in mathematics. In fact, it was clearly demonstrated that all the enabling and simplifying features 
of the CESE method are rooted in three key ideas, i.e., (i) conservation laws be enforced in their natural 
space-time unity forms; (ii) the simplest stencil be used; and (iii) a solver be built from a non-dissipative 
core scheme. 

Moreover, it has been shown that the a scheme possesses some rather intriguing properties such as: 
(i) its two independent mesh variables it" and (u $ )" separately satisfy two decoupled three-level leapfrog 
schemes Eqs. (2.23d) and (2.24d), and (ii) it shares with the leapfrog scheme the same amplification factors, 
even though the a scheme and the leapfrog scheme have completely different origins and structures. These 
properties are tributes to the power of beautiful and deep physical ideas underlying its development. 

In addition, we also establish rigorously the deep relation between STI (i.e., PI) invariance and non- 
dissipatedness. It represents a good example of universal applicability of a fundamental physics idea. 

Perhaps the biggest surprise that occurred during the CESE development is the realization that a two- 
level explicit solver cannot be STI invariant if, corresponding to one physical conservation law, it has only 
one mesh variable and one equation per mesh point. It implies strongly that a discrete system may behave 
fundamentally different from a continuous system. 
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c APPENDIX 

c 

implicit real*8 (a-h, o-z) 
parameter (nxd=2000) 

dimension u (nxd) , un (nxd) , ux(nxd), uxn (nxd) , vo (nxd) , v(nxd), 

* vn (nxd) , vxo (nxd) , vx(nxd), vxn (nxd) 

c 

c Program ''lfa.for'' — a combined dual ''a'' and leapfrog solver 

c for the pure advection equation. The computation is carried out 

c assuming a sine-wave initial data and periodic boundary 

c conditions. The sine-wave is of unit wavelength. It can be shown 

c that both exact and numerical solutions must be spatially 

c periodic with unit wavelength too. 

c 

c u, un, ux, uxn arraies associated with the $a$ scheme. 

c vo, v, vn, vxo, vx, vxn arraies associated with the leapfrog 

c scheme, 

c 

c ux, uxn, vxo, vx, and vxn represent the numerical analogues of 

c normalized spatial derivatives during time marching. However, 

c the output contains the non-normalized values, 

c 

c it = no. of time-marching steps, 

c k = no. of spatial intervals per unit length, 

c cn = Courant number, 

c a = the advection speed, 

c 

it = 1000 
k = 101 
cn = 0 . 5d0 
a = l.dO 
c 

dx = 1 . dO/df loat (k) 
dt = cn*dx/a 
hdx = dx/ 2 . dO 
t = dt*dfloat (it) 
kl = k + 1 
k2 = k + 2 

pi = 3.1415926535897932d0 
tp = pi*2.d0 
pdx = pi*dx 

cns = (l.dO - cn**2)/2.d0 
cnp = (l.dO + cn)/2.d0 
cnm = (l.dO - cn)/2.d0 
c 

open (unit=8 , f ile=' for008 ' ) 

write (8,10) it,k,dt 

write (8,20) a,cn 

write (8,30) t 

do 100 j = 1 , k2 

tpx = tp*df loat ( j-1 ) *dx 

u ( j ) = dsin (tpx) 

ux ( j ) = pdx*dcos (tpx) 

v(j) = u ( j) 

vx ( j ) = ux ( j ) 
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100 


continue 
do 400 i = 1, it 
do 200 j = 2 , kl 

un ( j ) = cnm*u(j+l) + cnp*u(j-l) + cns*(ux(j-l) - ux(j+l)) 
uxn ( j ) = ( u ( j + 1 ) - u(j-l))/2.d0 - cnp*ux(j + l) - cnm*ux(j-l) 

if (i.eq.l) goto 150 

vn ( j ) = vo ( j ) - cn* (v ( j+1) - v(j-l)) 

vxn ( j ) = vxo ( j ) - cn*(vx(j+l) - vx(j-l)) 

goto 200 
150 vn ( j ) = un ( j ) 

vxn ( j ) = uxn ( j ) 

200 continue 

do 250 j = 2 , kl 
vo ( j ) = v ( j ) 
vxo ( j ) = vx ( j ) 

250 continue 

do 300 j = 2 , kl 
u ( j ) = un ( j ) 
ux ( j ) = uxn ( j ) 
v ( j ) = vn ( j ) 
vx ( j ) = vxn ( j ) 

300 continue 

u (k2) = u (2) 
ux (k2 ) = ux (2 ) 
v(k2) = v (2 ) 
vx (k2 ) = vx (2 ) 
u (1) = u (kl ) 
ux ( 1 ) = ux (kl ) 

V (1) = v ( k 1 ) 

vx ( 1 ) = vx(kl) 

400 continue 

tpat = tp*a*t 
do 500 j = l,kl 
x = df loat ( j-1 ) *dx 
tpx = tp*x 
ue = dsin (tpx-tpat ) 
uxe = tp*dcos (tpx-tpat ) 
eu = u ( j ) - ue 
ux ( j ) = ux ( j ) /hdx 
eux = ux ( j ) - uxe 

ev = v ( j ) - ue 
vx ( j ) = vx ( j ) /hdx 
evx = vx ( j ) - uxe 

write ( 8 , 40) 
write (8,50) x,ue,uxe 
write (8,60) u(j), eu, ux ( j ) , eux 
write (8,70) v(j), ev, vx ( j ) , evx 
500 continue 

close (unit=8) 


10 

format ( 

' it = 

00 

■H 

= ' , i8, ' dt =' 

, gi4 .7) 


20 

format ( 

' a =' 

,gi4.7, ' 

CFL = ' , gl4 . 7 ) 



30 

format ( 

' t =' 

, gi4 .7) 




40 

format ( 

' **** 

********, 

************** 

*************** 

***********') 

50 

format ( 

' x =' 

,gi4.7, ' 

ue =' , gl4 . 7 , ' 

uxe =' , gl4 . 7 ) 


60 

format ( 

' u =' 

,gi4.7, ' 

eu =' , gl4 . 7 , ' 

ux =' , gl4 . 7 , ' 

eux =' , gl4 . 7 ) 

70 

format ( 

' V =' 

,gi4.7, ' 

ev =' , gl4 . 7 , ' 

vx =' , gl4 . 7, ' 

evx =' , gl4 . 7 ) 


stop 

end 
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Figure 1.— A surface element on the boundary 
S(V) of an arbitrary space- time volume V. 
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2(a).— The space-time mesh. 
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Figure 2.— The SEs, BCEs, and CCEs. 
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Figure 3.— FV(j,n+l). 
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Figure 5.— The basic stencils of the a scheme and their STI images. 
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